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Abstract 

Under the assumption that some part of the observed highest energy cosmic rays consists 
of protons originating from cosmological distances, we derive bounds on the associated flux of 
neutrinos generated by inelastic processes with the cosmic microwave background photons. 
We exploit two methods. First, a power-like injection spectrum is assumed. Then, a model- 
independent technique, based on the inversion of the observed proton flux, is presented. 
The inferred lower bound is quite robust. As expected, the upper bound depends on the 
unknown composition of the highest energy cosmic rays. Our results represent benchmarks 
for all ultrahigh energy neutrino telescopes. 
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1 Introduction 



Ultrahigh energy cosmic ray (UHECR) protons with energies above the Greisen-Zatsepin-Kuzmin 
(GZK) cutoff, E>Eqzk = 4 X 10^^ eV, interact inelastically with the photons of the cosmic 
microwave background (CMB) and produce pions. This process results in a significant energy 
loss and an attenuation length of about 50 Mpc. If the highest energy cosmic rays are protons and 
originate from distances beyond that scale, one expects a sharp drop in the observed spectrum at 
around -Egzk PIJ- This GZK phenomenon also predicts a guaranteed ultrahigh energy neutrino 
flux, since the produced pions finally mainly decay into neutrinos These neutrinos are called 
GZK or cosmogenic neutrinos. 

There are a number of estimates [2jl3] and upper bounds |Hl5|l6] on the cosmogenic neutrino fluxes, 
the most recent results being found in Refs. [UlSlinillllllll] • We summarized these predictions in 
Fig. ^ Surprisingly, one finds two orders of magnitude uncertainty due to the huge differences 
between the individual results. Present and future experiments need a clear picture of this 
phenomenon. Therefore, we present in this Letter a systematic quantitative analysis of the 
minimal and maximal expected cosmogenic neutrino fluxes. The only assumption we make is 
that some part of the observed highest energy cosmic rays consists of protons from cosmological 
distances. We include all cosmological and observational uncertainties into our calculation. 

The inferred lower bounds on the cosmogenic neutrino fluxes turn out to be quite robust. These 
bounds are of particular interest. First of all, they represent benchmarks for all neutrino telescopes 
and cosmic ray facilities designed to be sensitive in the ultrahigh energy region (for a recent 
review, see Ref. |IZj). Moreover, the lower bound on the flux can be turned into an upper bound 
on the neutrino nucleon cross-section, if no quasi-horizontal or deeply-penetrating air showers 
are observed ^Sl- This knowledge gives important information about a possible enhancement of 
the cross-sections in the multi-TeV centre-of-mass energy regime ^H], which is expected in many 
standard model (SM) like or beyond the SM scenarios. 

Our Letter is organized as follows. In Section |2l we summarize how the propagation of protons 
through the CMB and the associated production of neutrinos can be described by means of 
propagation functions. Section El presents two techniques to give upper and lower bounds on the 
cosmogenic neutrino fluxes. In Section |3 we conclude. 

2 Propagation functions 

Protons propagating from cosmological distances lose their energies by three basic processes. 
Above the GZK cutoff, the dominant particle physics process is scattering on the CMB through 
pion production. The produced pions decay, giving rise to the cosmogenic neutrinos. For proton 
energies between 10^^ eV and the GZK cutoff, the dominant particle physics process is scattering 
on the CMB through e^e~ pair production. The expansion of the universe redshifts all the 
travelling particles, which is particularly important for protons produced at large distances. 

The propagation of the proton toward the earth can be easily described [T^ l^ by one single 
function Pp\p{E; Ei,r), which tells the probability that a proton created at distance r with energy 
Ei is detected on earth as a proton with energy above E. This probability function was calculated 
in Ref. [25 for a large range of E, Ei and r. More generally, including other particle species, the 
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Figure 1: Various predictions for the cosmogenic neutrino fluxes. We assume full mixing between 
the flavours [12], and the presented results correspond to one flavour, Fj,^ + Fp^, i = e, fi, t. Lower 
(upper) short dashed line: flux from Ref. for redshift evolution parameters (cf. Eq. (0)) ^ = 2, 
^max = 2 (n = 4, Zinax = 4). Lowcr (upper) long dashed line: flux from Ref. [S], assuming a 
maximum energy of Ernax = 3 X 1020(21) gV for the ultrahigh energy cosmic protons (cf. Eq. (j3))). 
Lower (upper) dashed dotted line: flux from Ref. P^, assuming n = 3 (n = 4). Lower and upper 
dotted line: flux from Ref. jlO] . assuming (Emax, a, "n.) = (1x10^^ eV,1.5,3) and (3x10^2 eV,1.0, 3), 
respectively, where a denotes the proton injection spectral index (cf. Eq. (jU). Boxes show the 
upper bounds obtained by combining Fly's Eye [13] and Agasa [13] limits on deeply-penetrating 
showers in different energy bins [T^]. The dotted band labelled by Auger represents the expected 
sensitivity of the Pierre Auger Observatory to z/,- + z^t, corresponding to one event per year per 
energy decade (TB] . 



effects of the propagation can be described [221111123] by Pb\a{E; Ei,r) functions, which give the 
expected number of particles of type b above the threshold energy E if one particle of type a 
started at a distance r with energy Ei. 

In this Letter, we assume that the sources of the ultrahigh energy protons (nucleons) are isotropi- 
cally distributed and can be described by a comoving luminosity distribution Cp{r, E^) of protons 
injected with energy E^ at a distance r from earth, which gives the number of protons per unit of 
comoving volume, per unit of time, and per unit of energy. With the help of the above propaga- 
tion functions, one can calculate the differential flux of protons (6 = p) and cosmogenic neutrinos 
(6 = vi, vi) at earth. Their number N}, arriving at earth with energy E per units of energy, area 
(A), time it) and solid angle (fi), can be expressed as 



F,{E) ^ 



dE dA dt dfi 
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Note, that this formula can be easily generalized to arbitrary source luminosity distributions. In 
the following, we make the usual assumption [71IHllUl lTm iTT| that the r and Ei dependences of the 
source luminosity distribution factorize, Cp{r,Ei) = p{r) Jp{Ei), and that the redshift evolution 
of the sources can be parametrized by a simple power-law, 

Cp{r, Ei) = po (1 + Zi^r))"- 6{z - ^min) ^'(^^max " z) Jp{Ei) , (2) 

where the redshift z and the distance r are related by dz = (l + z) H{z) dr. We use the expression 
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to relate the Hubble expansion rate at redshift z with the present one. Uncertainties of the latter, 
EQ = h 100 km/s/Mpc, with h = (0.71 ± 0.07) 

^agi -"^^^ affect our results significantly. 

In Eq. Q, and ^a, with Qm + = 1, are the present matter and vacuum energy densities 
in terms of the critical density. As default values, we choose Qm = 0.3 and = 0.7, which is 
favoured today. Our results turn out to be quite insensitive to the precise values of the cosmolog- 
ical parameters within their uncertainties. The universe is homogeneous for distances above the 
GZK scale (50 Mpc). The events between 10^^ eV and the GZK cutoff are believed to come from 
cosmological distance. Therefore, it is physically motivated to use the ansatz Q for the evolution 
of the luminosity of sources in addition to the pure redshifting (n = 0). We set minimal and 
maximal redshift values by Zj^i^^ and Zmax- These parameters exclude the existence of nearby and 
early time sources. We take = 0.012, corresponding to r^in = 50 Mpc. The effects due to a 
change in ^Tmax can be compensated by a change in n. Therefore, we fix Zmax = 2 in the following 
and study the dependencies on n only. 

We neglect the effects of possible magnetic fields. They just deflect the trajectories of the protons 
and, thus, only increase the path length (synchrotron radiation of protons can be neglected). As 
long as the correlation length of the magnetic fields is smaller than the gyroradius of the protons 
(which holds for an anticipated magnetic field strength of ~ nG), this effect is on the percent 
level, smaller than other uncertainties. 

The details of our calculation of the Pf,^a{.E] Ei,r) functions for protons, neutrinos, charged lep- 
tons, and photons will be published elsewhere I2S]- In short (see also Ref. |E]), we calculated 
Pb\a{E; Ei, r) in two steps, i) First, the SOPHIA Monte-Carlo program [23] was used for the sim- 
ulation of photohadronic processes of protons with the CMB photons. For e~^e~ pair production, 
we used the continuous energy loss approximation, since the inelasticity is very small (~ 10"'^). 
We calculated the Pb\a functions for "infinitesimal" steps (1 -r- 10 kpc) as a function of the redshift 
z. ii) We multiplied the corresponding infinitesimal probabilities starting at a distance r{z) down 
to earth with z = 0. 

The determination of the propagation functions took approximately one day on an average per- 
sonal computer. The advantage of the formulation of the spectra in terms of the propagation 
functions is evident. The latter have to be determined only once and for all. Without the use of 
the propagation functions, one would have to perform a simulation for any variation of the source 
luminosity distribution Cp{r, Ei), which requires excessive computer power. Since the propagation 
functions are of universal usage, we decided to make the latest versions of —dPb\a/dE available 
for the public via the World- Wide- Web URL www.desy.de/~uhecr . 

3 Bounds on the cosmic neutrino flux 

In this Section, we present two techniques to derive robust upper and lower bounds on the 
cosmogenic neutrino fluxes. 

The first technique assumes an E'" power-like injection spectrum for the protons, with some 
maximal cutoff energy -Emax- A comparison of the spectrum after propagation with the obser- 
vations allows to determine the confidence region for the power a and for the redshift evolution 
index n. Since the propagation of the protons leads to neutrino production, we may then infer 
the neutrino fluxes suggested by the different (a,n) regions. 



4 



The second technique is based on an inversion of the observed proton flux around the GZK cutoff 
with the help of the proton's propagation function. Though the inverted proton spectrum has 
non-negligible uncertainties, the inferred fluxes of the cosmogenic neutrinos are rather stable. We 
study the sensitivity of the resulting flux on the cosmological evolution parameter n. 

The lower bounds for the neutrino fluxes obtained by these two techniques are in complete agree- 
ment. Since the post-GZK events can only be taken into account in the second method, the 
corresponding upper bound is larger than the one obtained by means of the first method. 

3.1 Power-like injection spectrum 

We assume an E^"" power-like injection spectrum for the protons, 

Jp{Ei) = Jo E^°' ^(-Emax — Ei) , (4) 

up to -Emax; the maximal energy which can be reached through astrophysical accelerating processes 
in a bottom-up scenario. The normalization factors Jq of the injection spectrum (0)) and po of the 
source distribution in Eq. (j2)) are fixed by the observed cosmic ray flux. The predicted differential 
proton flux at earth, Eq. (Q), with this injection spectrum, is compared with the observations. A 
fitting procedure gives the most probable values for -Emax, « and n. 

We quantify the goodness of our results by statistical methods. Our analysis is similar to that of 
Ref. [TT]. It was carried out in two basic steps. 

i) First we determined the number of experimentally observed events in a given energy bin by 
converting the published values of the cosmic ray flux. This had to be done, since the UHECR 
collaborations give their results for the observed flux in a binned form, whereas the number 
of events in a given bin is integer and follows the Poisson distribution. We analyzed the results 
from different experimental settings separately and performed the analysis for the two most recent 
results from the AGASA j^H] and HiRes |27j collaborations. In the low energy region, there are no 
published results available from AGASA and only low statistics results from HiRes-2. Therefore, 
we included the results of the predecessor collaborations - Akeno [2H] and Fly's Eye 29^ - into 
the analysis. With a small normalization correction, it was possible to continuously connect the 
AGASA data with the Akeno ones and the HiRes-1 monocular data with the Fly's Eye stereo 
ones, respectively (cf. Fig. El (left)). The normalization was matched at E = 10^^'^ eV for both 
cases. 

ii) We determined the 2-sigma confldence regions in the a-n plane for each -Emax- In order to 
do that, we checked the compatibility of different (a,n) pairs at a given -Emax with the observed 
data. In this analysis, we used the energy range between E^ = IQ^'^-'^^^^-^ eV and E^ = 10^" eV. 
The data in the bins above 10^*^ eV are not used, since events beyond the GZK cutoff can not 
be explained by just one power law injection spectrum with spectral index a. This obvious 
statement can be formulated quantitatively. Using z^i^ = 0.012 (which reflects the observation 
that there are no UHECR sources within rmin ~ 50 Mpc), we flnd that the data above 10^° eV 
are incompatible on the 3-sigma level with a pure power law flt in the energy region between 
E_ = 1017-2^18.5 gv and E+ > lO^o eV (see also Ref. jSOl). 

The compatibility of a given [a, n) pair with the observational data was checked as follows. For 
some speciflc (a,n) pair, the expected number of events in individual bins is calculated (A = 
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Figure 2: Two-sigma confidence regions in the power-law index ot and the redshift evolution index 
n plane. Results obtained with the Akeno + AGASA data are shown by a solid curve, whereas 
the results from Fly's Eye + HiRes data are given by the dotted line. The triangle (square) 
represents the best fit values for AGASA (HiRes). The analysis used the data between energy 
bins E_ = 10^^-2 eV and = lO^o eV. Other parameters were E^ax = 3 x lO^^ eV, z^in = 0.012, 



{Ai, Xr}, where the Aj's are non-negative, usually non-integer numbers, and r is the number of 
the bins). The probability distribution in the i-th bin is given by the Poisson distribution with 
mean Aj. The r dimensional probability distribution P(k) is just the product of the individual 
Poisson distributions (here k = {ki, ...kr} is a set of non-negative integer numbers). It is easy to 
include also the ~ 30% overall uncertainty in the energy measurement of the experiments into 
the -P(k) probability. According to the r dimensional probability distribution, the experimental 
result s = {si,...Sr} (where the Sj's are non-negative, integer numbers), has a definite, though 
usually very small probability P{s). The (a,n) pair is compatible with the experimental results 
at the 2-sigma level if 

^ P(k)<0.95. (5) 

k|P(k)>P(s) 

The best fit is found by minimizing the sum on the left hand side. This technique is equivalent 
to the technique for a large class of problems. 

For E_ = 10^^-2 eV and E+ = lO^o eV, our best fit values are ^^ax = 3 x lO^^ eV, a = 2.57, 
n = 3.30, for AGASA, and E^ax = 3 x lO^i eV, a = 2.50, n = 3.80, for HiRes^ Figure H displays 
the 2-sigma confidence regions in the a — n plane with -Emax = 3 x 10^^ eV for both experiments. 
Figure ini (left) shows our best fits to the Akeno + AGASA (top) and to the Flys's Eye -|- HiRes 
(bottom) UHECR data. Fig. El (right) shows the resulting neutrino fluxes per flavour, assuming 
full mixing at arrival at earth [1^ . 

For larger -Emax, these confidence regions are unchanged. If we lower -Emax, they decrease and 
disappear at -Emax = 3 x 10^° eV for AGASA and -Emax = 1 x 10^° eV for HiRes, respectively. In 
order to get bounds on the cosmogenic neutrino flux, we vary -Emax? and n, within their 2-sigma 
allowed values for both experiments. As there is still no consensus about the origin of UHECRs 
in the range from 10^'' eV to ~ 10^^ eV, we perform the analyses using both = lO^'''^ eV 
and E- = 10^^'^ eV, respectively. The only difference, in the latter case, is that n is no more 

^For similar analyses, with Zmin — 0, see Ref. |3J 
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Figure 3: Left panel: Ultrahigh energy cosmic ray data with their statistical errors (top: com- 
bination of Akeno and AGASA data; bottom: combination of Fly's Eye and HiRes data). The 
best fits between E_ = lO^'''^ eV and = 10^° eV, using the power-hke injection spectrum, are 
given by the solid lines. The 2-sigma variations corresponding to the minimal (dotted) and max- 
imal (dashed) fluxes are also shown. Other parameters of the analysis were -Emax = 3 x 10^^ eV, 

■^min = 0.012, and ^max = 2. 

Right panel: Neutrino fluxes per flavour, F^^ + Fp^, i = e,^,T. The "best" predictions for the 
neutrino spectra are given by the solid lines. The 2-sigma variations corresponding to the minimal 
(dotted) and maximal (dashed) fluxes are also shown. 
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Figure 4: Lower bounds on the neutrino fluxes, F^^ + Fp^, £ = e, fi, r, using a power-hke proton 
injection spectrum, starting the fit at E_ = 10^'^'^ eV (dotted) or E_ = 10^^-^ eV (dashed). In the 
latter case, n is not constrained and we used n = 3. Both are compatible with the lower bound 
obtained by the propagation inversion technique using n = 3 (solid). 
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constrained. As we will see in the next subsection, for a reasonable upper bound we have to 
assume that also the post-GZK events are protons from cosmological distances. This is clearly 
beyond the possibilities of a power-like injection spectrum. Therefore we conclude this subsection 
with a presentation of the inferred lower bound in Fig. 01 

3.2 Propagation inversion 

The basic strategy of this subsection can be summarized as follows. Protons from cosmological 
distances contribute to the observed high energy cosmic rays (e.g. above ~ 10^^ eV). By means 
of the propagation function of the proton, one can determine the proton injection spectrum. This 
spectrum contains protons above the GZK cutoff, which results in cosmogenic neutrinos. Using 
a "minimal" or "maximal" proton contribution to the high energy cosmic rays, we can derive 
model-independent lower and upper bounds on the ultrahigh energy neutrino fluxes. 

Usually, the propagation function is not invertible. A given detected spectrum Fh{E) can be 
produced by several source luminosity distributions Cp{r,Ei). As an illustration for this non- 
invert ibility, one may think of two (unphysical) extreme cases. It could be that the whole detected 
spectrum is produced by nearby sources and no energy loss takes place. In this case the injection 
spectrum is the same as the observed one. In the other extreme case all the protons are produced at 
some large redshift -Zmax- They lose their energy when they propagate through the universe and the 
injection spectrum contains much more high energy events than the observed one. Nevertheless, 
an unambiguous inversion can be carried out by fixing the r dependence by some physical choice. 
We shall again assume an isotropic source luminosity distribution of the form characterized 
by a redshift evolution parameter n and -Zmin/max- 

From Eqs. ((TJ) and we find 

F,{E)= dE,G,\p{E,E,)Jp{E,), (6) 

where Gb\p{E,Ei) is the space integral of Eq. ((H). Instead of calculating the integral over Ei, 
usually one performs a summation over the energy bins, 

F,{E) = J2AE,G^p{E,E,)Jp{Ei), or F, = G^pJp, (7) 

where the second equation is the short-hand notation for the matrix-vector multiplication. With 
the help of Gf,\p, it is straightforward to invert the observed proton spectrum and to determine 
the resulting cosmogenic neutrino spectrum (in this case the particle type of "6" is "neutrino" ) , 

F,. = Gu\p Gp|p Fp , (8) 

where Fp is the vector notation of the observed proton flux Fp{E). 

The inversion procedure has an additional complication due to the fact that the observed spectrum 
in the high energy region has only a few events. The lack of large statistics results in significant 
statistical uncertainties. In order to take into account these effects, we used a Monte-Carlo 
and, for both experiments, we generated 10^ hypothetical observed spectra, compatible with the 
experimental results. Applying Eq. (jH)) to these spectra, one obtains the most probable cosmogenic 
neutrino flux, together with its statistical uncertainties. 
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Figure 5: Bounds on the cosmogenic neutrino flux per flavour. Lower bounds use the fact that the 
accumulation in the normalized proton spectrum is a consequence of the GZK mechanism. Upper 
bounds assume that the whole spectrum, even the part above the GZK cutoff, originates from 
uniformly distributed proton sources with isotropic luminosity distribution. The region with the 
dashed boundary represents the maximal and minimal neutrino fluxes for which the cosmological 
evolution parameter was allowed to change between and 6. The region with the solid boundary 
is obtained by a cosmological evolution parameter n = 3. Boxes show the upper bounds obtained 
by combining Fly's Eye [121 and Agasa ^1] limits on deeply penetrating showers in different 
energy bins [TSj. According to these limits, strong cosmological evolution {n > 3) for the sources 
of the post-GZK events is already excluded. 



Inverting these spectra, one obtains the most probable injection spectrum with its statistical 
uncertainties. Propagating these spectra through the universe, one infers a neutrino spectrum 
with its statistical uncertainties. 

In order to enhance the different features of the observed spectrum, one usually multiplies it with 
the energy to some power, E"^, where 7 = 2-^3. In this "renormalized" spectrum, one observes 
an accumulation of events just around the GZK cutoff, as can be seen in Fig. |3] (left). A possible 
physical origin of this effect is apparent. Protons above the GZK cutoff lose their energy quite fast; 
however, as soon as they reach an energy around the GZK cutoff, their energy loss degrades very 
much. Thus, their relative number is larger. Other particles would produce such an enhancement 
at other energies and protons from nearby distances (< 50 Mpc) would not produce such an 
enhancement at all. The apparent observation of this accumulation suggests, therefore, that the 
observed spectrum around the GZK cutoff is dominated by protons and that their sources are at 
cosmological distances. In our ^^minimar scenario, only this part of the spectrum is assumed to 
be given by protons from cosmological distances. Since there are very good indications that this 
part of the spectrum is really a result of the GZK processes, the neutrinos, produced in the same 
processes, provide a robust lower bound on the cosmogenic neutrino flux. As a lower bound on 
the guaranteed cosmogenic part of the neutrino flux, it also serves as a robust lower bound on 
the total ultrahigh energy neutrino flux. 

There is no conventional astrophysical explanation for the observed comic rays events beyond the 
GZK cutoff. Even their particle composition is unknown. In our "mazzmaf scenario, we assume 
that they are all protons created isotropically in the universe. Thus, their injection spectrum was 
large enough to survive the GZK cutoff and produce the observed spectrum. Through the GZK 
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mechanism, also these protons produce neutrinos. Since in this scenario we assume that all of 
the observed particles are protons (and calculate the associated cosmogenic neutrino flux), there 
is no more room for additional cosmogenic neutrinos. The bound we obtain is a robust upper 
bound on the cosmogenic neutrino flux (however, not necessarily an upper bound on the total 
neutrino flux). 

In our inversion procedure, we used the observed cosmic ray spectrum between E_ and E+. We 
changed between 10^^'° eV and 10^^'^ eV. Our results are rather insensitive to this choice. 

In our ^^minimar scenario, we obtained the lower bound on the neutrino flux. In order to check 
the sensitivity of our result on E^, we lowered E^ from 10^" eV to 10^^'^ eV. The change of the 
neutrino fluxes is marginal for Ej^ values between 10^° eV and 10^^'^ eV. Taking an unphysically 
small value as low as 10^^'^ eV, the change in the spectrum is smaller than the uncertainty of the 
expected sensitivity of the Auger experiment at Ej^ = 10^^ eV (on the "renormalized" flux figure, 
the experiment is most sensitive at this energy, cf. Fig. The GZK process suppresses the 
spectrum extremely effectively above E ^ 10^'' eV. Deviations from an expected GZK suppressed 
spectrum starts to be statistically significant above E f» 10^° eV. Therefore, we used E^ = 10^° eV 
in our analysis, similarly to our analysis based on a power-like proton injection flux. The lower 
bound for n = 3 agrees well with the one obtained with the other technique (cf. Fig. H}. 

In our ^^maximar scenario, we used the whole observed spectrum up to the highest observed 
AGASA event in the E = 10^°'^ eV energy-bin. In principle, there could be other neutrino 
sources than the GZK process, therefore the total neutrino flux might be even larger. However, 
since our maximal cosmogenic neutrino flux almost reaches the experimental limits (cf. Fig. El, 
there is not much room for additional neutrinos. Note, that along with the cosmogenic neutrino 
flux there is also a cosmogenic photon flux from 7r° decay. After propagation, the photon flux 
associated with our upper bound in Fig. may be in conflict with the EGRET observation of the 
diffuse gamma ray flux [S2[IS1- For a detailed analysis, one would need the photon propagation 
functions ^B] as well. 

Our results on the neutrino flux bounds are summarized on Fig. El We used n = 3 cosmological 
evolution. The variation due to the change of this parameter between and 6 is also shown. 
Present [12] and expected (TE] experimental bounds by AGASA, Fly's Eye and Auger, respectively, 
are also presented. 

4 Conclusions 

In this Letter, we presented a systematic quantitative analysis of the minimal and maximal 
expected cosmogenic neutrino fluxes. The only assumption we made was that some part of the 
highest energy cosmic rays is due to protons from cosmological distances. We used two techniques. 
One of them assumed an power-like injection spectrum for the protons with some maximal 
cutoff energy -Emax- We compared the predicted spectrum after propagation with the observed 
one and determined the neutrino spectrum which was produced during the propagation. The 
second technique was based on an inversion of the observed proton flux around the GZK cutoff. 
The prediction for the lower bound on the cosmogenic neutrino fluxes were rather stable. 

These lower bounds are of particular interest. They represent benchmarks for all neutrino tele- 
scopes and cosmic ray facilities designed to be sensitive in the ultrahigh energy region. 
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